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Abstract 

Several variants are presented of a linear-in-parameters least squares formulation for determin- 
ing the transfer function of a stable linear system at specified frequencies given a finite set of 
Fourier series coefficients calculated from transient nonstationary input-output data. The basis 
of the technique is Shinbrot’s classical method of moment functionals using complex Fourier 
based modulating functions to convert a differential equation model on a finite time interval 
into an algebraic equation which depends linearly on frequency-related parameters. 
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1. INTRODUCTION 

Methods for determining the transfer function of a stable linear system from input-output 
data include correlation and spectral analyses, as well as the direct sinusoidal measurements. 
Each of these “nonparametric” identification techniques require either a statistical stationarity 
assumption on the data, or a periodic steady state condition to be established, before initiating 
calculations of the transfer function at pertinent frequencies. Excellent summaries of these 
methods, as well as the analysis of noise effects and finite data lengths, can be found in 
Astrom [1], Ljung [2], Soderstrom and Stoica [3], and Unbehauen and Rao [4]. Notwith- 
standing noise considerations, long data lengths may be required in order to achieve good 
accuracy due to the stationarity or steady state assumption. By contrast, a method is proposed 
here that utilizes the frequency content in short data lengths in order to set up a least squares 
estimation of the transfer function at selected frequencies. Since short data lengths are used 
there is no assumption of steady state operation or stationarity of the data, though there must 
be present sufficient energy content in the data at the specified frequencies in order to avoid 
nondegeneracy in the least square estimate. The basis of the technique is the classical Shin- 
brot [5] method of moment functionals, also known as the modulating function technique, 
using complex Fourier based modulating functions. A forerunner of this approach can be 
found in Pearson and Lee [6] which utilized real valued Fourier based modulating functions, 
i.e., commensurable sinusoids, for the parameter estimation of linear differential systems. 
Therein also can be found a discussion of the background of this method with a listing of 
references. However, this paper appears to represent the first use of modulating functions in 
the context of the “nonparametric” system identification problem. Several variants of a 
deterministic least squares estimation of frequency-related parameters that underlie the transfer 
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function will be developed below. 

2. Least Squares Formulations 

Consider a stable linear system with input u (t ) and output y(t) modeled on a finite time 
interval by the differential operator equation: 

A(p)y(t)=B(p)u(t)+e(t) (1) 

where (A(p)Ji(p )) are polynomials in the differential operator p=d/dt of degree less than or 

equal to an a priori integer n , and e (t ) represents the effect of modeling errors. The problem 
considered here is to estimate the transfer function G (J to )=B (j t o)M (j to) at a finite set of fre- 
quencies {£co 0 , * =1 ’ 2 ' ' Af }, where 0) 0 is a user selected “resolving” frequency and M a 
chosen integer, given the input-output data \u it ),y (t )] over a finite set of time intervals 
{[tiJi+T], i=l, ■ ■ N J. 1 These time intervals are each chosen of length 7’=2jt/co 0 and need not 
necessarily be disjoint. However, a certain degree of independence in the data collected over 
the different [0,7] time intervals is necessary in order to avoid degeneracies in the least 
squares estimate to be discussed below. Understandably these degeneracies are more likely to 
be avoided in normal operating records if the intervals are disjoint. In addition to the upper 
bound on the system order n , the DC value of the transfer function is assumed given or can 
be measured from the step response, i.e., G (0)=B (0)M (0) is presumed known a priori. If 
this is not the case, then the estimated transfer function can be scaled by the parameter G (0). 

The Shinbrot method of moment functionals is a technique for converting a differential 
equation on a finite time interval into an algebraic equation by the use of “modulating func- 

1 If the system bandwidth (0g is known, then choosing (M ,to 0 ) such that M tt) 0 =co 5 will cover 

the bandwidth at the knots k(. 0q, k= 1,2 • • M. 
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tions”. As introduced by Shinbrot [5], <J >(f) is a modulating function of order n relative to a 
finite time interval 0<f <7 if it is sufficiently smooth and satisfies the end point conditions: 

<^>(0 )=4> (,) (7)=0, 1=0,1 • • n-l (2) 

where 4> (,) (r)=p‘<t>(r). Clearly, modulating functions can be constructed in many different 

ways. 2 Here a specific set of complex valued Fourier based modulating functions is defined in 

a way that is conducive to solving the problem at hand, viz., let 

$ m (t)=e jm< * / (l-e j<s *) n , 0<t<T =2 k/(D 0 (3a) 

m = 0,1 • • M 

define a set of modulating functions of order n with respect to the time interval [0,7]. 
Equivalently by the binomial expansion, each such function is representable by 

<UO=£^ (m+ * )<U °'. (3b) 

k=0 

where b k is defined in relation to the binomial coefficient by 



The first representation (3a) makes evident the fact that each <]> m (f) is indeed a modulat- 
ing function of order n , i.e., (2) is satisfied, 3 while the second representation (3b) implies that 
calculating linear functionals defined by each <t> m (r) on a set of functions specified over [0,7] 
will entail calculating the Fourier series coefficients of these functions at the frequencies fcco 0 > 
k=m/n+l • • m+n, m=0,l • • M. In turn, these coefficients can be calculated efficiently by 
DFT/FFT methods which provides an important motivating factor for this analysis. This will 

2 See discussion in Pearson and Lee [6]. 

3 Notice that any modulating function of a fixed order is automatically a modulating function of 
any lower order relative to the same time interval. This property facilitates the formulation for any 
system of order less than or equal to the upper bound n . 


be discussed further below. The important property of the functions defined in (3) is con- 


tained in the following 4 

Modulation Property. Let Pip) be a differential operator of order at most n, i.e., a 
polynomial in p=d/dt of degree <n , and z (r ) any sufficiently smooth function defined 
on [0 ,T). Then the modulation of P(p)z(t) with <|> m (r) over [0 ,T) satisfies 

T m+n 

fo m (t)P(p)z(t )dt = X b k-m p Hk Q>0 )Z-k (5 ) 

0 k=m 

where Z k is the k th harmonic Fourier series coefficient of z(t), i.e., 

T 

Z t =Jz(r)< (6) 

0 

Note that owing to the end point constraints (2) satisfied by each <|) m (0 function, none of the 
boundary point derivatives z^(0) or z^‘\T) appear in (5). This is crucial to the ensuing 
analysis and, in fact, represents a primary reason for employing the modulating function tech- 
nique. 

2.1. Formulation 1 

A direct application of the above property to the problem posed involves rewriting the 
differential operator model (1) in the equation error form followed by the modulation with 

<|> m (r); thus, 5 


4 Proof of this property is given in the Appendix in order to proceed directly with the develop- 
ment. 

5 The process of going from the model (1) to equation (7) can be viewed as a projection from a 
space of functions on [0,7] down into a finite dimensional space spanned by the modulating functions. 
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fo m (r) \a (p )y (t )-B (p)u(t)]dt =fo m (t )e (t )dt. (7) 

o L o 

In view of (5) the preceding equation is equivalent to 

m+n r "1 «+", 

x h-m U W* C0 0 )y_, -B (-jk a> 0 )t/_* I = Ysh-mE-k- (8) 

k=m L J k=m 

Define the real and imaginary parts of the polynomials (A (jk <a 0 )Jl (jk w 0 )) as follows: 

A (jk co 0 )=a k +j (3^ , B (jk co 0 )=y* +j 8* (9) 

and collect these together to form the 4x1 “parameter” vector: 

a* 

e* = £ • (>°) 

A. 

Also, define as follows the 2x4 data matrix \|/*(0 in terms of the real and imaginary parts of 
the jfc^* harmonic Fourier series coefficients of the input-output data corresponding to the time 
interval [r, ,r, +T ] , / =1,2 • • N : 

Y c k (i) Y k (i ) -U c k (i) -Ui(i) 

V*0)= _ Y s (i) yc (l) ^ (0 -{/£(;) 

The notation for the entries in (11) is explained by 

r t 

Y k (i ) = jy (f +r,- )cosfc to 0 r dt , Y s k (i ) = Jy (r+r, )sin£ co 0 r dt (12) 

o o 

and similarly for (U£(i),U k (i)). Then the real and imaginary parts of the equation error (8) 
can be collected into the following real valued 2x1 vector equation which serves as the start- 
ing point for a least squares estimation: 

m+n 

2>*-mV*(0©* =£m0') 03) 

k=m 



m=0,\ ■ ■ M , i=l,2- -N. 
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The equation error vector in (13) is related to the Fourier coefficients of the original equation 
error by 

m+n 

O' ) — S k-m 

k=m 

Values of the transfer function G (jk co 0 )=5 (Jk co 0 )/A (jk (Oq) are seen from (9) and (10) to 
be related to the parameter vector 0* by the real and imaginary part relations: 


m) 

W).‘ 


ReG (jk o) 0 )= 


a *+P* 


ImG (jk co 0 )= 


<**8*-P*V k 


or equivalently by the magnitude-phase relations: 


(14) 


“"*■ (15) 

Starting from a presumed knowledge of the DC value G (0), which implies that 
0(p[/l (0),0^ (0),0f is known, equation (13) can be rewritten in the standard regression equa- 
tion format to estimate the parameters Q k , £=1,2 • • M+n given the data over a sufficient 
number of [r i ,r ( +7'] ) intervals, i=l,2 • • N. A consideration of this equation reveals the fol- 
lowing: 

1) The frequency range covered by the parameters in (13) is (M+n)c o 0 . Hence, if it 
is desired that the transfer function estimate cover a frequency range about 25% 
greater than the system bandwidth co fi at a resolution cOq, a choice in (M ,co 0 ) such that 


(M+n)co 0 = 1.25 g) b (16) 

reflects this objective. 

2) Counting unknowns in (10) and (13), the total number is 4 (M+n). Since each 
equation in (13) is of dimension 2, counting equations suggests that the total number 
N of [ti,ti+T] intervals should satisfy: 2(M +l)N>4(M +n), i.e., N>2(M+n )/(M+\). 
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However, the equations in (13) are partially decoupled with respect to the index m. 
Therefore, it seems best to estimate (0j, • • 0 n ) at the first stage, which corresponds to 
m - o. This means that there are 4 n unknowns for the first stage requiring that N 
satisfy N>2n. Thereafter, the number of unknowns is just 4 for each stage 
corresponding to m= 1,2 • • M , which implies N > 2 assuming that the preceding esti- 
mates are used in each succeeding stage. This kind of “bootstrapping” of the least 
squares estimation facilitates keeping the number of unknowns to a modest level at 
each stage. 

3) The two row vectors comprising y*(0 in (11) are seen to be mutually orthogonal 
for each k and i suggesting a maximal degree of independence for these equations in 
utilizing the information content in the data. This is a direct result of the Fourier 
nature of the underlying formulation. 

Discussion: The above development shows that it is possible to formulate a linear-in-the 
parameters least squares estimation problem for parameters (10) that underlie (via (14) or 
(15)) the transfer function G(jk to 0 ) at each k th harmonic frequency. The input-output data 
can be time-limited and transient, but must have sufficient energy content at the specified fre- 
quencies to avoid degeneracies in the least squares solution. Apart from being highly non- 
linear, the relations (14) and (15) involve the difference between parameter related quantities, 
e.g., a k y k -$ k b k , whose values may be large for large k. This aspect of the formulation por- 
tends a potential source of error magnification which is alleviated by the formulation of the 


next section. 
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2.2. Formulation 2 

Given that [u(r),y(r)] satisfies the model (1) on a [0,7] time interval, it follows that 
[«(r)j(t)] also satisfies the model 

A (-p )A ip )y (t )=A (-p )B (p )u (r )+A {-p )e{t). (17) 

Choosing a set of modulating functions of order In to accommodate the upper bound on the 

highest degree differential operator in (17) and modulating this equation with the m th member 

of this set, the following projected equation error results which is analogous to (8): 

m+2n r I m+2n_ 

E h-m U (jku 0 )A Hk(O 0 )Y_ k -A (jku 0 )B (-jk(a 0 )U_ k 1= £ b k _ m A (jku 0 )E_ k (18) 

k=m L J k=m 

where b k is defined by, cf. (4): 

h = <-»‘ ft") 

Noting that 4 (jka 0 )A (-jk 0 i Q ) is real while A (Jk (£> Q )B (-jk (0 0 ) is complex, define real quanti- 
ties (a^.otfc.P*) by the relations: 

a k =A (jk co 0 M (~jk O) 0 ), a* +j P/t = A (Jk “o ) 5 (rfc w o) ( 19 ) 

and collect these into the 3x1 parameter vector 0* defined by 

^k 

0 * = (x k . ( 20 ) 

p* 

Also, define the 2x3 data matrix \y k ( i ) by 

Y c k (i) -U k (i) Utii) 

V * 0)= Y s k (i) -U s k (i) -U c k (i) 

where the notation for the entries in (21) is the same as defined in (12). Then the real and 

imaginary parts of the projected equation error (18) can be represented by the following real 

2x1 vector equation: 
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m+2n 

Z Vm'MOO* =e m (0 (22) 

k=m 

m= 0,1 -M, i= 1,2 - N. 

Equation (22) can be rewritten into the standard regression equation format for setting up 
the least squares estimate of the parameters (0j, • • 0^+2 «) based on the data over [r, ,t[+T], 
i=l,2 • • N. Again, presumed knowledge of the DC value implies that 0 O =4 (0)[A (0),S (0),0]' 
is known or, if not, the resulting transfer function estimate can be scaled by the parameter 
G (0). Here the estimates of the transfer function are related to the parameters by the real and 
imaginary part equations (as found from (19) and (20)): 

ReG (jk C 0 q)= — — , ImG (jk co 0 )=-— (23) 

a k a k 

or equivalently by the magnitude-phase relations: 

|G (jk co 0 )p = — -ij— ' , /g (Jk 03 q) = -tan -1 ^-. (24) 

a k a * 

Consideration of the least squares formulation in this case leads to the following: 

1) The frequency range covered by the parameters in (22) is (M+2n )co 0 ; hence, the 
guideline (analogous to (16)) for choosing the pair (M,0) 0 ) in this case is 

(M+2n )a>o = 1.25c0fl . (25) 

2) Counting unknowns in (20) and (22), the total number is 3 (M +2 n ) which would 

imply that the total number N of [t h ti+T] time intervals should satisfy: 
2A(A/+l)>3(M+2n). However, the partially decoupled nature of the equation (22) 
with respect to the m index suggests bootstrapping the solution from the first stage. 
Thus, for the initial stage (m= 0), N needs to satisfy N>3n. In the succeeding stages, 
N needs to satisfy N> 2 (since there are 3 unknowns and 2 N equations at each such 
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stage); this assumes that the preceding estimates are used at each succeeding stage 

(m= 1,2 • • M). 

Discussion: Comparing (23) and (24) with (14) and (15) reveals that the second formulation 
avoids the potential error magnification problem of differencing large estimated quantities in 
calculating the transfer function at high frequencies. However, the second formulation 
requires estimating 6 n unknowns at the first stage, i.e., the m=0 stage, verses 4n unknowns 
for the first stage of the first formulation. 


2.3. Formulation 2-Dual 

The dual to the formulation of the preceding section is to observe that a given pair 
\u (t ),y (t )] satisfying the model (1) on a [0,7’] time interval also implies that it satisfies (cf. 
(17)) 


B (-p )A (p )y (r )=B (- p )B(p)u(t )+B {-p )e (r ). (26) 

Again choosing a set of modulating functions of order 2 n, a development similar to that of 

the previous section leads to the real 2x1 vector equation 


m+2n 

52 ^k-m Vfc 0 — O' ) 

k-m 

m =0,1 * * M y 1=1,2 * * N 
where the data matrix \\f k (/ ) is defined by 


¥*( 0 = 


U£(i) - Y c k (i ) Y£(i) 
U£(i) -Y s k (i) -Y k (i) 


(27) 


(28) 


and the real 3x1 parameter vector 0* is defined by 
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with the entries in (29) defined by 



b k 

a* 

-P* 


( 29 ) 


b k =B (jk co 0 )B (- jk co 0 ) , +j $ k = A (Jk to q)B ( —jk co 0 ) . (30) 

Relations between the transfer function and the parameters in this case are found to be 

<Xl blr “Pit b k 

ReG (Jk co 0 )= , ImG (jk co 0 )= — (3 1 ) 

a* 2 +p * 2 «*+P* 

for the real and imaginary parts, or for the magnitude-phase: 

2 a 

P (jk co 0 f = bk ~ - , / 3 (jk co Q ) = -tan” 1 -^ . (32) 

Comparing (21) and (28) verifies the duality of the two formulations by virtue of the 
interchange of input and output. Note that each formulation has the same total number of 
unknowns - in general. However, the dual formulation has the potential advantage of reduc- 
ing the total number of unknowns in the event of a priori information on a lower degree 
numerator polynomial than denominator polynomial in the transfer function. For example, an 
“all pole” model means that b k =(B (0)) 2 is known for all k, i.e., a total of 2 (M +n ) unknowns 
verses 3(M+2n) for the previous case. The formulation leading to (27) can easily be 
modified to reflect this consideration. 


3. Conclusions 

Three formulations of a linear-in-the parameters least squares estimation have been 
presented for determining the transfer function of a linear system at specified frequencies 
given transient nonstationary input-output data. While some comparisons have been noted in 
the discussions following each formulation, a clear indication of the pitfalls and advantages of 
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each will have to await a thorough simulation study including the effects of noise. Such a 
study is currently underway. 
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5. Appendix 

To verify the modulation property (5) it is sufficient to consider the differential operator 
P(p)=p i for a fixed i<n. The result for a general n th degree polynomial Pip) then follows 
by superposition. Thus, for sufficiently smooth z(f) on [0,7’] and § m (t) defined in (3), the 
left side of (5) in this case is 

T T 

[4 ) m (t)p i z(t)dt = fz it)p l $ m (t)dt (33) 

o o 

where integration-by-parts has been used i times taking into account the boundary conditions 
(2) possessed by each <f» m (r) function. Substituting the representation (3b) into the right side 
of (33) and changing the index of summation verifies (5) for P(p)=p‘ as purported. 


